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Abstract 



We investigate in some detail an alternative simulation strategy for 
lattice field theory based on the so-called worm algorithm introduced by 
Prokof'ev and Svistunov in 2001. It amounts to stochastically simulating 
the strong coupling expansion rather than the usual configuration sum. A 
detailed error analysis and an important generalization of the method are 
exemplified here in the simple Ising model. It allows for estimates of the 
two point function where in spite of exponential decay the signal to noise 
ratio does not degrade at large separation. Critical slowing down is practi- 
cally absent. In the outlook some thoughts on the general applicability of 
the method are offered. 



* e-mail: uwolff@physik. hu-berlin . de 



1 Introduction 



Since about a decade [T] there has been a development in the condensed matter 
community where a novel approach to the simulation of statistical models was 
pioneered. One replaces the original partition function given by a sum over certain 
field configurations by its untruncated series representation, usually the strong 
coupling or hopping parameter expansion. For any finite but arbitrarily sized 
system this is regarded as an alternative statistical system that is simulated by 
a novel Monte Carlo technique. Originally the method has emerged in quantum 
models that are 'made classical' by using the Trotter formula which leads to world- 
line type formalisms. In [2] the approach was used for classical spin systems 
in two and three dimensional- In our opinion this idea is very useful also for 
investigating systems of interest to lattice field theory and elementary particle 
physics and it might be applicable in suitably generalized form even to fermions 
and gauge models. Ultimately it could perhaps enable steps toward the long sought 
generalization of [3] beyond infinite coupling^. 

In this publication that is planned to be the first of a series, we concentrate 
our numerical experiments on the 2D Ising model as a prototype system. While 
the extension to abelian scalar bosonic theories in any D seems straight-forward - 
some cases can already be found in [2] other generalizations will be non-trivial. 
The purpose of this paper, apart from some reformulation, is on the one hand to 
introduce an important generalization of biasing the sampling of the series. On 
the other hand a more detailed investigation of the behavior of some observables 
of interest and their comparison with standard methods is made. It should have 
become clear above that we are not discussing just a new Monte Carlo algorithm 
but an exact reformulation where physical observables have completely different 
estimators, variances and autocorrelations which deserve study. Another study 
of the worm algorithm in the Ising model was recently presented in [3]. There a 
detailed study of the Monte Carlo dynamics is given refining results in [2]. While 
we confirm the practical absence of critical slowing down on our lattices, we here 
do not attempt to determine dynamical critical exponents. 

For a first closer look at the method we define 



Here x, y, u, v denote sites of a cubic periodic lattice in arbitrary dimension D. The 
outer sum is over Ising configurations {s(x) = ±1}, in the exponent we sum over 

1 I would like to thank Urs Wenger for drawing my attention to this paper. 

2 Similar hopes have been recently expressed by S. Chandrasekharan in his talk at Lattice 
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all links /. This is the numerator of the two point function that may be written as 



G(x-y) = (s(x)s(y)} = ^\, (2) 

Zj [Z, Z) 

with an arbitrary lattice site z. Generalizing [2] we consider a new partition 
function 

Z=^2p-\v-u)Z(u,v). (3) 

u,v 

The pairs of sites u, v are now the 'phase space' and < p{x) < oo is a weight 
function demanded to possess the lattice periodicity Constant positive factors in 
p will be irrelevant and we adopt the normalization condition at the origin 

p(0) = 1. (4) 

Now the two point function is given by 

G(z) = p(z)%^, (5) 

where the double bracket ((. . .)) denotes averages with respect to (jHJ) 



for observables v). Note that any dependence of mean values on the weight 
p is canceled in (GO), but statistical errors will depend on it. 

To actually simulate (j3J) we need Z(tt, v) as weights which are not available in 
closed form. Inserting (fj]) and enlarging the phase space to {u, v,s} would lead 
to oscillating over-all weights and would not be useful. Instead here the strong 
coupling expansion may be inserted. In a somewhat sketchy notation, that will 
become more precise later, we put 

Z(u,v) = k (tanh/3) d(7) . (7) 

7Sr('U,u) 

Here the strong coupling expansion of the Ising model is written in powers of tanh (5 
and k is an unimportant constant. The set of strong coupling graphs contributing 
to ([T]) is named T(u,v) and for each element 7 of the set 0^(7) denotes the degree 
in tanh (3 of this contribution. Of course T(u, v) includes disconnected graphs that 
may wind around the torus etc. If we insert this representation into ([3]) we may 
consider {it, v, 7} as phase space points that may now be sampled probabilistically 
to estimate ((. . .)). The very important result of [2] is that it is almost trivial 
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to design an ergodic update from joint essentially local moves of {u, v,^y} under 
which observables as those appearing in (jSJ) are even (almost) free of critical slowing 
down. Local here means in particular that in each step 7 is deformed only in a 
minimal way close u or v. A detailed description of the moves will follow. For the 
graphs of the usual partition function ('vacuum' graphs) alone an ergodic update 
procedure would be more complicated, it is easier in the enlarged setup. At the 
same time (jSJ) shows that the extended moves produce information on the two 
point function if we record the occurring v — u in a histogram. 

It may seem surprising that local moves can (almost) eliminate critical slowing 
down. But this prejudice is of course based on the idea of sampling spin configura- 
tions that are long distance correlated. Here as one approaches criticality, bigger 
strong coupling graphs will be important and apparently the algorithm 'knows' 
this and may even be guided in addition by choosing p appropriately. 

There also is an important difference with regard to internal symmetries, the 
spin flip Z(2) in the Ising case. Spin-flip odd observables cannot be addressed 
anymore. The contributing graphs reflect the symmetry manifestly and the field 
variables that transform under the symmetry have been summed over. This should 
not be seen as a drawback: In the ideal situation of knowing the exact solution 
this would be the same. 

The plan of the paper is as follows. In 2. we derive the dimer form of the 
Ising model that explicitly parametrizes all strong coupling graphs. In addition 
the update for this system is defined. In 3. some numerical results are collected 
and discussed. We end in 4. on conclusions and an outlook, where we offer some 
speculations on the possibilities to generalize this reformulation and its simulation 
in various directions. 

2 Ising model demonstration 

2.1 Dimer formulation 

For a single Ising bond the trivial identity 

e Ps(x)s(y) = CQsh p [tanh(P) s(x)s(y)] k (8) 

fc=0,l 

holds. Using it on each link with independent dimer or bond variables {ki = 0, 1} 
we rewrite ([T|) (up to the factor 2~ Nx inserted for convenience) as 

Z(u,v) = (cosh/3)^2-^ nt tanh ^) II s(x)} kl *(u)s(v) (9) 

{k,s} I x€dl 
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with the number of links (sites) iVj (N x ). The notation I = (xy) has been replaced 
by the boundary set dl = {x,y}. 

Next, the spins are summed over and leave behind the constraint 

2 ~ Nx E IM s(x)]^s(u)s(v) = Q(k; it, v) G {0, 1}. (10) 

{s} I x&dl 

The result B factorizes into local constraints with 

«(*;») = { J L c '*»*' ) = even (ii) 

and their complements 

0(*;jO = 1-0(A;;I/) (12) 

combining to 



e(k;u,v) = ( n ^v))><{l^%ell,v) =V else. < 13 > 



In words: If u = v holds, the number of dimers touching at any site has to be even. 
For u v , these two sites have to be surrounded by an odd number of dimers with 
all other sites being even. 

The previous steps coincide with those made in a duality transformation [5]. 
There, usually with no insertion (or u — v), one would solve the constraint by 
expressing fc; by (the dual of) the discrete 'exterior derivative' [0] of a D — 2 'form' 
on the dual lattice, yielding self-duality for D = 2, the dual gauge theory for 
D = 3, etc. For the algorithm [2], however, we stay with the constrained {hi}, in 
a sense intermediate between the original Ising model and its dual. 

The partition function of the enlarged ensemble now reads 

Z= Y e fo^)e-MS.». (14) 
^ p(v-u) 
u,v,{h} ^ v ' 

where the dimer sum together with the constraint uniquely label the graphs 7 G 
T(u, v) and the total dimer number equals 0^(7). The coupling has been traded for 
the dimer chemical potential /1, 

tanh/3 = e- M . (15) 

Mean values of observables A(k; u, v), which may now also depend on k h are given 
by 

<W>=4 E A(k;u,v)^^le-^. (16) 
Z ^-^ p(v — u) 

u,v,{ki} 



5 



Simple diagnostic observables for first Monte Carlo experiments are the internal 
energy, or equivalently the average nearest neighbor correlation 



T,J2 G ^ ( 17 ) 

and the susceptibility 



E D 

V 



X 



where we sum over all directions v and v are the corresponding unit lattice vectors. 
The transcription of \ is obvious, 



X- 1 = N f U ' V \ v (19) 

Thus, in particular for p = 1, the fraction of coinciding u = v has a simple 
interpretation and vanishes for (5 ^ j3 c in the thermodynamic limit. 
The corresponding translation for E is 

1 {{p{v - u)6\ u - v \,t)) 
E = 2D ((M) • (20) 

Alternatively it may be related to the mean dimer occupancy in the subset u = v 
by differentiating hi Z(z, z) with respect to (3 



J>u,v N Yll k , , 

E = e~» + 2 sinh // " ' " . (21) 

((*«,«)) 

2.2 Prokof 'ev-Svistunov worm algorithm 

With the algorithrr|§ of [2] we sample the ensemble (I14p . The name of the algorithm 
presumably derives from the fact that the constraint Q(k; u, v) ^ requires a line 
of active dimers (ki = 1) connecting the ends u, v of the worm. However the 
connecting sequence of links is not unique, the worm is 'fuzzy' so to speak with 
fixed head and tail only. We therefore prefer the name PS-algorithm honoring its 
inventors. 

In any case the PS algorithm exploits the fact that we can base an ergodic 
Monte Carlo algorithm for the ensemble flHj) on merely the following two types of 
simple elementary steps that we apply in alternating order: 



'Originally with p = 1 only. 
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• type I: we pick one of the 2D nearest neighbors of v with equal probability 
and call it v' and the connecting link I. The proposed move v — > v' with 
the simultaneous adjustment k\ — > 1 — h\ is accepted with the Metropolis 
probability 



otherwise the old configuration is maintained. The practical implementation 
is largely based on precomputed tables. It is not necessary to also move u 
due to translation invariance. 

• type II: if we encounter a configuration u, v, {k{\ with u = v we 'kick' u = v 
together to a randomly chosen other lattice site with unchanged {ki} with 
the probability < pkkk < 1- For u ^ v we do nothing in this step, which is 
the dominant case. 

Each such I-II step requires 0(1) operations, independent of the lattice size. We 
call this compound a micro-step. For a sequence of 2Ni micro-steps ergodicity can 
be shown. Any two configurations may be connected with non-zero probability by 
the 'worm dismantling' all active dimers of the first configuration and then 're- 
building' the second one. The jumps allow to move from one connected component 
to the next. As usual in ergodicity proofs, this is in general of course a rather 
academic process! Following [2] we employ pkkk = 1/2. Some brief numerical 
experiments with values 0.3 and 0.7 showed that this has hardly any influence. 
Metropolis acceptance rates in steps I were reasonably far from the extreme values 
and 1 in all our simulations below, usually between 0.5 and 0.6. In pE] it has 
been observed that a correct and efficient algorithm can even be based on steps of 
type I alone (pkick = 0) which we have however not tried here. 

3 Numerical demonstrations 

We group together N x micro-steps calling this an iteration, during which we ac- 
cumulate a sub-histogram of occurring separations v — u. This may be interpreted 
as measuring the whole set of 1-bit observables {S V - U)Z } - one for each separation 
z - after each micro-step, mostly implicit zeros, of course. Thus we 'always mea- 
sure' and never give away any information. Such an iteration has a computational 
complexity comparable to a sweep in standard algorithms. We separately record 
the contributions to primary observables from each such iteration. These records 
are treated as a usual time series and an autocorrelation analysis as described in 
[7j is made. 

Most physically interesting observables are of the derived type, nonlinear func- 
tions, typically ratios, of primaries, see (0). Details on their error estimation, 




(22) 
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including the definition of integrated autocorrelation times for such functions, can 
also be found in [7] . We only repeat here that in our definition \/2rbt is the ratio 
between the true and the 'naive' error. Although this is wide-spread, there are 
alternative definitions in the literature that may coincide with this one only for 
large not relevant here. 

The tests below are in D = 2 using lattices of size T x L, but they are expected 
to be representative for larger D, too. 

3.1 Simulations at f3 c with p = 1 

Results for x m the form f)19p and E from both (1201) and (|2T|) are reported. For 
each lattice size we have performed 10 6 iterations after equilibration^. The values 
are compatible for instance with [8] and our results with integrated autocorrelation 
times are listed in tabled! We confirm the practically complete absence of critical 
slowing down for the observables under study. 



T = L 


L 7 'VX 


'int,x 


£(J2D} 


r int,£;j2a 


£(EQ) 


T int,i?l2U 


16 


0.9166(28) 


0.813(14) 


0.7263(10) 


0.567(5) 


0.7267(5) 


1.94(4) 


32 


0.9153(33) 


0.815(14) 


0.7171(9) 


0.535(6) 


0.7171(3) 


2.13(4) 


64 


0.9117(38) 


0.860(15) 


0.7116(8) 


0.517(5) 


0.7120(2) 


2.41(5) 


128 


0.9196(44) 


0.928(17) 


0.7096(7) 


0.508(3) 


0.7097(1) 


2.60(6) 


256 


0.9102(49) 


0.933(17) 


0.7089(7) 


0.504(3) 


0.7083(1) 


2.81(7) 



Table 1: Values and integrated autocorrelation times for simulations at the critical 
point (3 = ln(l + y/2)/2 and various lattice sizes. 



It is instructive to compare the achieved precision with that of the single cluster 
simulations in combination with standard spin estimators that were reported in 
[8]. The energy E is a useful monitor for simulations but non-universal and not 
really of physical interest. We hence concentrate on x which contains long distance 
contributions. We first compare the growth of T int x here with the one quoted in 
[8] for the single cluster algorithm (1C). While both algorithms are very fast, 
the PS simulations seem to have even smaller and less rising Ti n t iX . Does this 
automatically mean that PS is preferable for measuring %l As we estimate in a 
completely different ensemble also the variance can be quite different. In [H] errors 
are quoted for x together with the number of single cluster steps performed. Using 
the average cluster size (that actually equals x) h is an eas Y exercise to predict 
the errors for a single cluster run of 10 6 steps per spin, a CPU effort similar to the 
one in this study. Taking L = 256 as an example, the result is that the error on x 
from a single cluster simulation would be about 6 times smaller. 

4 Hcre and below we generously spend 10% of the total statistics for equilibration. 
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Due to the simplicity of the estimator ( fl9l) we immediately know the variance 

v x -i,rs = «(*«,„ - X' 1 ) 2 )) = X- 1 ~ X- 2 « X" 1 . (23) 

Given the known final errors and the fact that we effectively measure after each 
micro-step, we can also conclude the integrated autocorrelation time Ti ntiX in time 
units of micro-steps. It ranges from about T intiX ~ 9 for L = 16 to T intiX 53 for 
L = 256. This means that an effectively independent estimate for x _1 arises after 
moving dimers (modifying the strong coupling graphs) in only a small part of the 
volume. The point is that we are discussing an integrated autocorrelation time and 
not the slowest mode. Our observable couples for instance only weekly to remote 
parts of the graph. Nevertheless, as discussed above, taking into account also 
the variance, no miracle in terms of efficiency has happened! In [I] quite similar 
observations have been made for x- in their comparison with the Swendsen-Wang 
cluster algorithm PS fares however somewhat better. 

We end this subsection by visualizing a typical contributing graph in figure [U 




Figure 1: A typical contribution for T = L = 64, (5 = (3 C without bias, p = l. The 
blobs are at u and v, lines show k\ — 1. 
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3.2 Biased PS: making use of p 



One use of the generalization of arbitrarily biasing the ensemble (I14p with p ^ 1 
may be seen as follows. We assume to be able to guess an approximate form 
of the two point function (up to a factor) and use this for p. Then §5$) shows 
that the Monte Carlo produces an approximately ^-independent correction factor 
that turns our guess into the exact correlation. This means, up to fluctuations 
and imperfections of our guess, all our counting-bins of occurring v — u will have 
approximately equal filling. Apart from autocorrelation effects - which will be 
found to be minor - this means that the relative statistical error of correlations 
measured in this way will be c/ yN with the number of iterations N and essentially 
the same constant c for all separations. This implies a constant signal to noise ratio 
independent of the size of the signal, a rather unusual situation. 

We first tested this mechanism in the critical model simulated before. The 
actual method to set p was: In a first simulation with p = 1 where we spent 
about 10% of the planned total statistics, we have measured G after some equi- 
libration. Although not mandatory, we symmetrized the observed histogram 'by 
hand' over the geometric lattice symmetries: reflections along the two directions 
and, if T — L holds, cubic rotations. These symmetries were, of course, only 
violated by the statistical fluctuations. With this p we then performed the main 
simulation^. We have indeed observed almost flat ((5 V - U:Z )) and compatibility with 
the known results. Nevertheless there is no significant gain in precision which was 
not unexpected: At criticality, in a finite volume, the two point function does not 
decay exponentially but only by relatively moderate factors. The real strength 
of the possibility of biasing toward larger separations is expected in the massive 
disordered phase (3 < (3 C which we study next. 



3.3 Mass gap in the disordered phase 

In many applications one is interested in low-lying eigenvalues of the transfer 
matrix. We pick a Euclidean time direction and resolve the lattice site label into 
two-vectors x = (xq,xi) labeling sites on a T x L torus and work in lattice units, 
a = 1 (integer x u ). We want to study the temporal decay of correlations at fixed 
spatial momentum pi 

G Pl (x ) = J2G(xo,x 1 )e- i ^ x \ (24) 



It is not difficult to see that in the dimer representation it is given by 

((p(t,v 1 -u 1 )5 V0 _ U0 ^ -^-^)) 

«<w> 



g pi (t) = ' 1 ^gp^f a (25) 



In more difficult cases one could think of a sequence of improving approximations. 



10 



In the present study we restrict ourselves to observing the mass gap and probe 
with pi = only. Moreover, we found by a short test that for this observable 
there is no disadvantage in letting p only depend on time. Then the correlation 
simplifies to 

So(t) = p(t) %^ . (26) 

We now need a guess of the time-slice correlation and use this for p and have the 
Monte Carlo again only work out the corrections. It is sufficient to histogram only 
the time separations and in addition count coincidences u — v, if the normalization 
is needed. Due to symmetry we may combine t and T — t. 

Below, instead of showing the correlation itself, we shall consider the time 
dependent effective mass. Taking care of time periodicity we define m(t + 1/2) by 
numerically solving for m in 

goft + 1) cosh(m(T/2-t-l)) <07\ 

— ^ / ; — = n — v^r, r; — > m > 0, (J ^ t < 1 / 2 (27) 

G {t) cosh(m(T/2 - £)) ' ; K J 

which for mt 3> 1 is expected to plateau at the asymptotic mass-gap of the transfer 
matrix. This finite volume mass gap has a weak L-dependence and approaches 
the true gap for large mL. Here we are not interested in this second limit but just 
work around mL « 5 and below (3 C . 

For a test-bed with L = 64, T = 256, (3 = 0.42 we show in figure [2] the effective 
mass from several simulations. The horizontal lines are the exactly known result^ 
from [9] 

m cx [L = 64,/3 = 0.42] = 0.0841370.... (28) 

Each of the plots corresponds to 10 6 steps per lattice site. The uppermost panel 
shows the standard situation in a favorable case. We have carried out 10 6 sweeps 
of the Swendsen-Wang algorithm and have measured the standard spin-estimator. 
Due to our large statistics and a well isolated spectrum, we see a clear plateau 
forming before the errors explode at larger separation. As usual in such estima- 
tions, while the signal Go vanishes exponentially the variance remains essentially 
constant. This could in principle be improved by cluster estimators, but the signal 
to noise ratio would still vanish exponentially with t, only at a slower rate. 

In the next panel from above we see the result of a PS simulation with p = 1. 
The general impression is similar as before. The errors are somewhat larger at 
small t and smaller at large t. For the next simulation we have biased for a flat 
histogram using an analytical guess p oc cosh(m(T/2 — t)) with a mass close to the 



6 I would like to thank Martin Hasenbusch for advise with the Ising literature and for a C-code 
that evaluates the finite volume mass-gap. 
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Figure 2: Plots of the effective mass (times L) versus separation t in correla- 
tion length units. From top to bottom: standard spin simulation, p = l,p — 
cosh(m(T/2 — t)), p = cosh((m + 6/T)(T/2 — t)). Errorbars are one sigma. 
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expected one. We are rewarded with a completely uniform error independent of t. 
In additional experiments we tried to over-sample long distances and have used 

poccosh[(m + 2n/T)(T/2-t)] (29) 

which exponentially favors v o — uq = T/2 over vq = uq by a factor exp(n). The 
bottom plot in figure [2] shows the case n — 3. We indeed see even smaller errors 
at large t (about a factor 0.6) where systematic errors are smallest, a reversal of 
the standard situation where one usually has to compromise between systematics 
and statistics. The three PS runs took equal CPU time and the cluster simulation 
in our implementation was a factor two longer, all on a single PC. 

For extremely non-uniform choices of p we expect to eventually run into a prob- 
lem, because u, v have to travel over the whole lattice including all separations for 
ergodicity. If one approaches such a case gradually, as we did varying n, this should 
not set in abruptly but autocorrelation times should start to rise. For the case 
above, and even for n > 3, we have hardly seen such effects except that sometimes 
the summation window [7j for the autocorrelation function, in particular for E, 
had to be increased to capture its tail. All observed integrated autocorrelation 
times of effective masses were very close to 1/2 in units of iterations. In addition 
we notice that in the PS data the errorbars jump around the exact line while the 
spin simulation results exhibit longer wavy structures. This indicates that the PS 
data from different time separations are statistically less dependent. Averaging 
the effective mass over a range or performing a fit in some suitable window would 
probably reduce the error significantly. If we take the deviations from the known 
constant for say mt > 2 and the determined errors to form a x 2_va l ue ) we find 
reasonable values close to the number of data points entering. We do however not 
attempt to optimize the mass gap extraction any further here in this demonstration 
of principles. 

4 Conclusions and Outlook. 

We have investigated the worm algorithm of Prokof 'ev and Svistunov in the case 
of the tanh/3 expansion of the Ising model. We emphasize that this involves a 
reformulation of the problem as a different statistical sum rather than just a new 
simulation algorithm to produce thermalized spin configurations. We confirm [2] 
and |4J to the effect that critical slowing down in practice is no issue for this 
method, at least in the case studied here. The important generalization of biased 
sampling (p ^ 1) allowed us to estimate the two point function in the disordered 
massive phase with a signal to noise ratio that does not decay with distance. Even 
the contrary can be achieved. 
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We expect that these features immediately generalize to other scalar models 
with abelian field variables, e.g. 4 , Z(iV) clock models, U(l) XY-model, Potts 
models in arbitrary dimension. A high precision study using the PS algorithm 
for the phenomenologically important XY model in D = 3 was already presented 
in [TO]. In [TT] the elimination of the sign-problem for bosons in a non- vanishing 
chemical potential by the PS-method is demonstrated. For correlations beyond 
two point functions one could consider more insertion points, perhaps with par- 
tially restricted domains. For non-abelian spin models like 0(N), RP(iV), CP(iV), 
SU(iV) x SU(iV) the expansion and constraints are clearly more involved but hope- 
fully still manageable. It should perhaps be kept in mind that for the continuum 
limit one is not confined to the standard actions of these models. An immediate 
generalization will be the world-line or loop-gas formulation of lattice fermions as 
presented in [12J. A similar representation can be written for fermions in more than 
two dimensions. It will be all decisive for the future usefulness of the approach 
whether the sign-problem can be handled there. In two dimensions it is essentially 
absent and successful worm simulations have in fact already been reported [T3] . 
This will allow for the simulation of the O(N) invariant Gross- Neveu model where 
along the lines of p2] we have encountered some problems with the measurement 
of interesting correlation functions. Another important generalization would be 
lattice gauge theory. There the insertion of spin fields obviously should be gen- 
eralized to Wilson loop insertions and gauge invariance will be manifest term by 
term. In this way it will hopefully be possible to efficiently generate the strong 
coupling graphs. The closed surfaces of the vacuum graphs will be embedded in a 
larger class of graphs with surfaces with boundaries as discussed to some degree 
in [IT] . Also here the abelian case will be much simpler than for example SU(2). 
First efforts towards a nonabelian dual formulation (without so far employing a 
PS-type algorithm) can be found in [H]. The ultimate dream, of course, would be 
an efficient alternative simulation technique for full QCD. 

Acknowledgments: I would like to thank Willi Rath and Urs Wenger for dis- 
cussions and Oliver Bar in addition also for useful comments on the manuscript. 
Support by the Deutsche Forschungsgemeinschaft (DFG) in the framework of SFB 
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